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Abstract. The stories told in this paper are dealing with the solution of finite, infinite, and 
biinfinite Toeplitz-type systems. A crucial role plays the off-diagonal decay behavior of Toeplitz 
matrices and their inverses. Classical results of Gelfand et al. on commutative Banach algebras yield 
a general characterization of this decay behavior. We then derive estimates for the approximate 
solution of (bi)infinite Toeplitz systems by the finite section method, showing that the approximation 
rate depends only on the decay of the entries of the Toeplitz matrix and its condition number. 
Furthermore, we give error estimates for the solution of doubly infinite convolution systems by 
finite circulant systems. Finally, some quantitative results on the construction of preconditioners 
via circulant embedding are derived, which allow to provide a theoretical explanation for numerical 
observations made by some researchers in connection with deconvolution problems. 
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0. Introduction. Toeplitz-type equations arise in many applications in mathe- 
matics, signal processing, communications engineering, and statistics. The excellent 
surveys ||, [l7| describe a number of applications and contain a vast list of references. 
The stories told in this paper are dealing with the (approximate) solution of biinfinite, 
infinite, and finite hermitian positive definite Toeplitz-type systems. We pay special 
attention to Toeplitz-type systems with certain decay properties in the sense that the 
entries of the matrix enjoy a certain decay rate off the diagonal. In many theoreti- 
cal and practical problems this decay is of exponential or polynomial type. Toeplitz 
equations arising in image deblurring are one example (since often the point spread 
function has exponential decay - or even stronger - compact support) [fj"9| . Kernels 
of integral equations also frequently show fast decay, leading to Toeplitz systems in- 
heriting this property (see e.g. Other examples include Weyl-Heisenberg frames 
with exponentially or polynomially decaying window functions |25f| (yielding biinfi- 
nite block- Toeplitz systems with the same behavior when computing the so-called 
dual window) , as well as channel estimation problems in digital communications pjj ] . 

Let C[_i i) be the set of all 1-periodic, continuous, real- valued functions defined 
on [—5, 5). For all / G C[_i i), let 



2 

a k = J f(uj)e 27!luk duj, fc = 0,±l,±2,. 
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be the Fourier coefficients of /. Since / is real-valued, we have = 5— fc. 

A Laurent operator or multiplication operator associated with its defining function 
f can be represented by the doubly infinite hermitian matrix L = [Lki\^'i = _ ao with 
entries Lu = ctk-i f° r k,l € Z. For all n > 1 let L n = [{L n )kiYkiL- n +i ^ e * ne Toeplitz 
matrix of size (2n— 1) x (2n— 1) with entries (L n )ki = afc_; for fc, Z = — n+1, . . . , n— 1. 
L„ is a finite section of the biinfinite Toeplitz matrix L. 

A Toeplitz operator with symbol f can be represented by the singly infinite matrix 
T = pfc;]?°,_ with Tfc; = afc_; for fe, Z = 0, 1, .... In this case we define T n = 
[(T n ) ki]k~[=o as the n X n matrix with entries (T n )kl — a-k-i for k, I = 0, . . . , n — 1. Of 
course L„ = T 2 „_i, but in what follows it will sometimes be convenient to use the 
notations L n and T n . 

As mentioned earlier, a crucial role throughout the paper plays the decay behav- 
ior of Toeplitz matrices and their inverses. Classical results of Gelfand et al. lead 
to a general characterization of this decay behavior for biinfinite Toeplitz matrices, 
see section |IJ Section ^| is concerned with the approximate solution of (bi)infinite 
Toeplitz systems using the finite section method. Explicit error estimates are derived, 
showing that the approximation rate depends only on the condition number of the 
matrix and its decay properties. In section || we analyze the approximate solution 
of convolution equations via circulant matrices. Finally, in section ^, we derive some 
quantitative results for preconditioning of Toeplitz matrices by circulant embedding. 
Among others, we provide a theoretical explanation of numerical observations made 
by Nagy et al. in connection with (non)banded Toeplitz systems. 

1. On the decay of inverses of Toeplitz-type matrices. It is helpful to 
review a few results on the decay of inverses of certain matrices. In what follows, if 
not otherwise mentioned, the 2-norm of a matrix or a vector will be denoted by ||.| 
without subscript. 

The following theorem about the decay of the inverse of a band matrix is due to 
Demko, Moss, and Smith Q. 

Theorem 1.1. Let A be a matrix acting on £ 2 (J), where J = {0, 1, . . . , N— 1}, Z, 
or N, and assume A to be hermitian positive definite and s -banded (i.e., A^i = if 
\k-l\ > s). Setn = IIAIIHA- 1 !!, q= ^1 and A = qi . Then 

\A^\<c\\ k -\ 

where 

c=p- 1 ||max{l,ii±^}. 

One notes that the inverse of a banded matrix is in general not banded, the 
type of decay changes when we switch from A to A -1 (although exponential decay is 
"almost as good" as handedness). This observation suggests to look at other classes 
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of matrices, for which the type of decay is preserved under inversion. This leads 
naturally to the following 

Definition 1.2. Let A = [Aj-^kjei be a matrix, where the index set isX = Z, N 
or {0, ...,N- 1}. 

(i) A belongs to the space £ 7j a if the coefficients Aki satisfy 

\A kl \ <ce-^ fc -^ / O r 7 ,A>0, 

and some constant c > 0. If A = 1 we simply write £ 7 . 

(ii) A belongs to the space Q s if the coefficients Am satisfy 

\A kl \<c(l + \k~l\y s fors>\, 

and some constant c > 0. 

The following result is due to Jaffard [ JL6| . 

Theorem 1.3. Let A : £ 2 (2) — > £ 2 (T) be an invertible matrix, whereX isl = Z,N 
or {0, ...,N- 1}. 

(a) If A £ £ 1 , then A^ 1 £ £ 7l for some 71 < 7. 

(b) If A £ Q s , then A' 1 £ Q s . 

For finite-dimensional matrices these results (and in particular the involved con- 
stants) should be interpreted as follows. Think of the n x n matrix A n as a finite 
section of an infinite-dimensional matrix A. If we increase the dimension of A n (and 
thus consequently the dimension of (j4 n ) _1 ) we can find uniform constants indepen- 
dent of n such the corresponding decay properties hold. This is of course not possible 
for arbitrary finite-dimensional invertible matrices. 



Theorem 1.3(a) shows that the entries of A -1 still have exponential decay, however 
A^ 1 is in general not in the same algebra as A, since we may have to use a smaller 
exponent. However in Theorem |l.3| (b) both, the matrix A and its inverse A -1 belong 
to the same algebra, the quality of decay does not change. 

From this point of view Theorem [01(b) is the most striking result. The proof of 



Theorem 1.3(b) is rather delicate and lengthy. For biinfinite Toeplitz-type matrices 
this result can be proven much shorter (and extended to other types of decay) by 
using classical results of Gelfand et. al. on certain commutative Banach algebras. The 
following theorem is a weighted version of Wiener's Lemma. It is implicitly contained 
in |1C| ], but since it may be of independent interest we state and prove it explicitly. 

Theorem 1.4. Let A — {a k i} be a hermitian positive definite biinfinite Toeplitz 
matrix with inverse A^ 1 = {a k {\- Let v(k) be a positive (weight) function with 

v(k + l) < v(k)v(l), 

such that 



00 

E 

k— — oc 



\a k \v(k) < 00. (1.1) 
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lim — = 1 and ]hn n y/v(n) = 1, (1.2) 

" T yu(-n) ™->°° 



n—>oo n. 



then 

oo 

\a k \v{k) < oo. (1.3) 



k— — oo 

In particular, 



if AeQ s for s > 1, then A- 1 e Q s ; (1.4) 
i/ A € £ 7) a /or < A < 1, iften A -1 e £ 7>A . (1.5) 



Proof. Since A is positive definite we have 

oo 

/H= ^ a fc e 2 ™^>0 (1.6) 

k——oo 

and by the properties of Laurent operators p"2| 

oo 

l//(w)= ]T a fe e 27rife -, where (A" 1 )*,, = a fc _ ; . 

fc — — OO 

We denote by W[u] the set of all formal series / = YlkL-oo a kX k for which 

oo 

11/11= M«(*)<°°- 

k— — oo 

It follows from Chapter 19.4 of [jlO| that W[v] is a Banach algebra with respect to the 
multiplication (discrete convolution) 

oo oo / oo \ 

fg= J2 *X l = J2 E a l-khjX l , 

1 — — OQ 1 — — OC \k= — oo / 

where / = J2k a kX k and g — J2k bkX k . By Theorem 2 on page 24 in |l(| an element 
of W[v] has an inverse in W[v] if it is not contained in a maximal ideal of W[v]. Any 
maximal ideal of W[v) consists of elements of the form (cf. Chapter 19.4 in fit! ) 

oo 
k— — oo 

where £ = pe 2vluJ with 



Pi < P < p2, 
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and 



Pi 



lim - 

n — *oo 7 



and P2 = lim 



Due to assumption (1.2) we get pi = P2 = 1, hence p = 1. Thus a necessary and 
sufficient condition for an element in W[v] to be not contained in a maximal ideal 



of W[v] is J2k a k 



^ for all lo. By assumption A is positive definite, hence 



f{uj) = Y^ k a k e 2 ™ kuJ > for all w and (|l_3|) follows. 

Statements ( |1.4| ) and ( |1.5| ) arc now clear, since in both cases we can easily find a 



weight function such that (1.1) and (|L2j) are satisfied. □ 

Remark 1.5. (i) Theorem 2.11 in || by Domar and Theorem V B in by Beurl- 



ing are closely related to Theorem 1.4. Their results are concerned with (non)quasi- 
analytic functions, for which they have to impose the more restrictive condition 



k=l 



\og[v(kx)] 
P 



< CO, 



for all 



on the weight function (called Beurling-Domar condition in p2[). For instance the 
function v(k) = exp( 1+ iog(|fc|) )' ^ ^ ^ satisfies condition (1.2), but not the Beurling- 
Domar condition. 

(ii) Using Theorem 8.1 on page 830 in p2| we can extend Theorem 1.4 to biinfinite 



block- Toeplitz matrices with finite-dimensional non- Toeplitz blocks (i.e., Laurent op- 
erators with matrix- valued symbol). These matrices play an important role in filter 
bank theory [^5j . 

(iii) Note that v{n) = exp(7n) does not satisfy condition (fb^), that is why we have 
to introduce an exponent 71 < 7 in order to estimate the decay of A -1 , cf. also The- 
lTy. However if A S £ 7 .a with A < 1, then condition 1.2 is satisfied and 



orem 



as 



we have seen - the decay of the entries of A 1 can be bounded by using the same 
parameters 7, A. 

2. Approximation of infinite-dimensional Toeplitz-type systems. Infi- 
nite Toeplitz systems arise for instance in the discretization of Wiener-Hopf integral 
equations or, more generally, in one-sided infinite convolution equations, see O]. Bi- 
infinite Toeplitz-type systems are encountered in doubly infinite (discrete) convolution 



equations, as well as e.g. in filter bank theory 25 or in the inverse heat problem |4j]. 
In order to solve these problems we have to introduce a finite-dimensional model. 

For let A : i 2 (Z) 1— > ^ 2 (Z) be a hermitian positive definite (hpd for short) biinfinite 
Toeplitz matrix given by {a k,l\W _oo- Let y — {j/fc}fcL_ oc G £ 2 (Z). We want to solve 
the system Ax — y. 

For n € N and y G £ 2 (Z) define the orthogonal projections P n by 



P n y = (. . . ,0, 0,y- n+1 , . . . ,y„_i,0,0, . . .). 



(2.1) 
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By identifying the image of P n with the 2n— 1-dimensional space C 2n 1 we can express 
the (2n — 1) x (2n — 1) matrix A n as 

An P n AP n , 

where we have used that P* = P. The n-th approximation af( n ) to a; is then given by 
the solution of the finite-dimensional system of equations 

A n x^ = 

where y( n > := P n y. 

If A is a singly infinite Toeplitz matrix and y £ £ 2 (N) , we proceed analogously by 
defining P n as 

P n y = (yo,J/i).--)J/n-i,0,0,...)< 

This approach to approximate the solution of Air = y is usually called the finite 
section method, cf. [ pd| . 

The first question that arises when considering this method is "does x^ converge 
to x?" . For the case when A is not hpd this question has lead to deep mathematical 
results. See the book and chapter 7 in Q for more details. For the case when A 
is hpd the answer is easy and always positive. To see this, recall that since A is hpd it 
follows that A n is also hermitian positive definite, see 
and ||(A n ) _1 || < for n = 1,2,.... Applying the Lemma of Kantorovich pllfl 

yields that (A„) _1 converges strongly to A~ x for n — ► oo, i.e., x^ converges to x in 
the ^ 2 -norm for any y £ € 2 (Z) (or for any y £ ^ 2 (N) if A is singly infinite). 

An important aspect for applications is if we can give an estimate on how fast 
x^ converges to x. It will be shown that the rate of approximation depends on the 
decay behavior and the condition number of the matrix. 

Theorem 2.1. Let Lx = y be given, where L = is a hermitian positive 

definite biinfinite Toeplitz matrix and denote x^ — L^y^ . 

(a) If there exist constants c,c' such that 

\a k \ < ce- 7|fc| and \y k \ < t/e -7 '*', 7 > (2.2) 

then there exists a 71 with < 71 < 7 and a constant c± depending only on 71 and on 
the condition number of L such that 

\\x-xW\\ < ce" 71 ". (2.3) 

(b) If there exist constants c,c' such that 

\a k \ < c(l + |fc|)- s and \y k \ < c'(l + \k\)~ s , s > 1, (2.4) 
then there exists a constant c\ depending only on the condition number of L such that 

\\x-xW\\ < cinC" 28 '/ 2 . (2.5) 



. Furthermore, \\A n \\ < \\A\\ 
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Proof. We have 



x™\\ = \\L~ x y - L-'y^W < W^My LL?y™\\ 

<\\L-i\\(\\y-y^\\ + \\(L n -L)L-^\\). (2.6) 



To prove statement (a) we note that by Theorem 1.3(a) there exists a 72 < 7 
such that {L~ l )t-i < C2e~ 72 l fc ~'l with a constant C2 depending only on 72 and on the 
condition number of L n . Since a(L n ) C [/ m j n , / max ] we get cond(L„) < cond(L) for all 
n. That means we can choose C2 independently of n. Write z^ n ' = (L n — L)L~ x y^ n > 
and note that 2! — for |fc| < n. Since the non-zero entries of (L n — L) decay 
exponentially, it is easy to show that there exists a 71 with < 71 < 72 such that 
ll z ^l < c 3 e~ 71 ™ for some constant c 3 . It is trivial that \\y — y^\\ also decays 
exponentially for n — ► 00. We absorb 1 1 1 1 j and the other constants in the constant 
Ci and get the desired result. 

For the proof of part (b) we proceed analogously to above by applying Theo- 
rem |1.3| (b) to conclude that 

|(i-V n) ) fe |<c 2 (l + fc)- s 

for some constant C2 depending only on cond(L) and on s. The norm \\y — y^ n '\\ can 
be estimated via 



\y-y (n) \\ 2 



00 „ 1 — 2s 

£ |y fc | 2 <2c£(l + fc)- 2s <2c J (l + x )-^dx<2c^—j, 

(2.7) 



|/c|>n k—n 



n-1 



similarly for ||z( n '|| where z^> :— [L n — V)L^y^ . Since all arising constants - 
absorbed in one constant c\ - depend only on cond(L) and on the exponent s, the 
proof is complete. □ 

Remark 2.2. Theorem [2.l| holds if we replace the system Lx = y by a singly 
infinite Toeplitz system Tx = y with corresponding decay conditions on T and y and 
approximate its solution by considering the finite system T n x^ n ' — y^ n \ 

The proof of Theorem |2.l| is essentially based on the fact that under appropriate 
decay conditions on L and y (resp. L n and jfW) L~ x and have similar decay 

properties. Thus, if one can show that L~ x has the same decay properties as 



one can use Theorem IA in order to generalize Theorem 2A to various other decay 
conditions. This may however not always lead to simple and closed-form expres- 
sions for || (L n — L)L~ 1 y'' n '\\, therefore I have restricted myself to the most frequently 
encountered decay properties. 



Example 1: We illustrate Theorem 2.1 by a numerical example. We consider 



Lx = y, where L is the biinfinite Toeplitz matrix with entries = (l+|fc|) s , k G Z for 
s = 2 and y consists of random entries having the same polynomial decay rate as the 
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actual error Hx-x'"' 11 
estimated error 




100 150 200 250 
n (size of circulant matrix is 2n+1 ) 



350 



Fig. 2.1. Actual approximation error and error estimate (2.5) from Theorem 2.1(b) for the 



system Lx = y, where L is a biinfinite hermitian Toeplitz matrix with polynomial decay. 



entries a&. To compare the error ||a:— x^\\ with the error estimate (2.5) we would need 
the true solution x. Since the solution of this biinfinite system cannot be computed 



analytically we compute the "true" solution of Lx = y by solving L r , 



r("o) 



y 



(no) fo r 



very large no (we choose no = 32768). Using ( |2.6| ) and ( |2.12| ) we can estimate that 
in the worst case \\x — a;( 32768 )|| w 10 -6 , so that a;( 32768 ) can mimick the true solution 
with sufficiently high accuracy for this experiment. 



Then we approximate this solution by the finite section method as in Theorem 2.1 



for n = 0, . . . ,350 and compute for each n the error \\x — x^ n '\\ as well as the error 



estimate in (2.5). Note that an explicit expression for the constant C\ in (2.5) is not 
known, we only know that it depends on the condition number of L. In this example 
we use ci = cond(L) (a different example may require a different choice). The result, 
illustrated in Figure [2.l| , shows that the asymptotic behavior of the error rate is well 
estimated by the given error bound. 

It is well-known that the product of two Laurent operators and the inverse of 
a Laurent operator (if it exists) is again a Laurent operator. This is of course not 
true for singly infinite or finite Toeplitz matrices (and this is one of the reasons 
which makes the "Toeplitz business" so interesting). Hence one may argue that the 
"canonical" finite-dimensional analogue of Laurent operators are not Toeplitz matrices 
but circulant matrices, since they also form an algebra. Thus for a given biinfinite 
hermitian Toeplitz matrix L with entries L^i — dk-l we define the hermitian circulant 
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matrix C n of size (2n — 1) x (2n — 1) by 



C n = 



ao ai ... a„_ 2 a n — i a n -i &n-2 
ai ap fli 



"l 



. . . ai 

« 2 

ai a 



(2.8) 



We also say that C„ is generated by {eife}£Z_„ +1 . 

Remark 2.3. C n does not have to be positive definite if L is positive definite, 
e.g. see ||. However - as pointed out in || - if L is at least in Wiener's algebra then 
one can always find an N such that C„ is invertible for all n > N. The faster the 
decay of the entries of L the smaller this N has to be. 

We can do even a little better and estimate how well the extrema of the defining 
function of L are approximated by the extreme eigenvalues of C n . 

Lemma 2.4. Let L be a biinfinite hermitian Toeplitz matrix with entries Lk,i = 



ttk-i where a = {a k } k x L_ 



id set f(uj) = J21 



ai-e 



2-Kiku 



Let C rn be the associated 



(2m — 1) x (2m — 1) circulant matrix with first row (ao, Si, . . . ,a m _i, a m -i, 
Denote the maximum and minimum eigenvalue resp. of C m by and 
(a) If L is n-banded with m > n, then 

nn 



AL m ax </ma X <A 1 t x + 2sin( 

A min - /) 



2(2m- l) 2 
7rn 



^-- 2 ™( 2 ,2 ra -l,= 



a||i, 



\a h. 



(b) If \a k \ < ce-^ forkeZ,c> 0, then 

2c 



\fn 



A (m) I < 

max I — 



l 



2 sin 



. / 7rm \ 
V2(2m- I) 2 J 



--ym 



e-i 

a similar estimate holds for |/ m i n — ^rnml- 

(c) If \a k \ < c(l + \k\)- s for k G Z,s > l,c> 0, then 

2c 



. ,ai). 

(2.9) 
(2.10) 

(2.11) 



A (r ' 



< 



1 



/ nm \ 
2 sin — --7 + 

V2(2m- 1 2 / 



2(2m 



(2.12) 



a similar estimate holds for \ / m i n — A 



(m). 



Proof, (a) : It is well-known that the eigenvalues of C m are given by 
,3*4*1/(2™-!) for/ = -m+l,..., 



m— 1 

E 

k— — m+1 



a k e 



1. 



(2.13) 



For case (a) that means they are regularly spaced samples /( 2 , ) of the function 



/, which immediately yields the left hand side of the inequalities (2.9) and (2.10). In 



order to prove the right hand side of (2^9) and ( 2.1C ) it is sufficient to estimate 

I 



max 



2m — 1 



2m - 1 
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where u> S [— 2 (2m-i) ' 2(2m-i) ] an< ^ ' = ~ m + 1, . - - , jtt, — 1. Define the sequence 
{°*}fcL_ m = a k if |fc| < n and a k = if |fc| > n. There holds 



max 



/(_!_)- /(-l±i^-) .max V 5fce 2-«/(2m-i) (e 2^/( 2m -i) _ 1} 



< max | e 2*i*W(2m-l) _]! V" | a I 
w,|fc|<n' z — ' 



k=-7 



< max 2|sin(7rfcw/(2m-l))|||a||i<2siii — ^ IMIi- (2.14) 

LO,\k\<n Z[Z']Jl, — 1 ) 



Relations and ( J2.10| ) follow now from this estimate. 



Statements (b) and (c) can be proved similarly by using 



/ -X- 

l/(^i)-*{ m) l< E klle 2 - ttw/(2ff, - 1) -i|+ E N- ( 2 - 15 ) 



k=-m+l |fe|> T 



and applying the corresponding decay properties to (2.15). □ 



Remark 2.5. The left part of inequality ( 2.10| ) reads A™ in > / m ; n . This im 



plies that Strang's preconditioner is always positive definite for n-banded hermitian 
Toeplitz matrices of size > (2n x 2n) with positive generating function. Hence the 
"sufficiently large n" -condition at the end of section 2 in |Q can be omitted. 

It is obvious that decay properties for circulant matrices cannot be defined in 
the same way as for non-circulant matrices. Hence, by stating that C m has, say, 
exponentially decaying entries, we mean that the generating sequence {a k }™~^ rn+1 
satisfies \a k \ < ce _7 l fc l, in which case the entries of C m will decay exponentially off 
the corners of the matrix (instead of off the diagonal) . In analogy to the theorems in 
section [l] it is natural to ask if the inverse of C m also inherits these decay properties. 
The following theorem shows that at least for m sufficiently large this is the case. The 
entries of C 1 " 1 uniformly approximate the entries of L^ 1 with an error rate depends 
on the decay properties and the condition number of L. 

Theorem 2.6. Let L be a hermitian positive definite Laurent operator with 
entries L k i = ojfe-j and let C m be the associated circulant (2m— 1) x (2m— 1) matrix 



as defined in (2.8). Denote the entries of L^ 1 by (L^ 1 )^ = {ak-i} and let (C m ) _1 
(if it exists) be generated by {f3 k } 1 ^~_ m , l ■ 

(i) If L £ £ 7 , then for sufficiently large m 

\a k -p k \<ce'^ m , (2.16) 

with < 71 < 7 and some constant c depending on cond(L) and 71. 

(ii) If A G Q s for s > 1, then for sufficiently large m 

\a k -/3 k \<c (2.17) 

s — 1 

with some constant c depending on cond(L). 
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Proof, (i): Set f(x) = EZ-oo a ke 2mkx and f m (x) = ET=- m+ i a k e 2mkx . By the 
properties of Laurent operators |12] {ak}kei, is given by 



1/2 



e 2 ™ kx dx, keZ, 



-1/2 



fix) 



(2.18) 



By Remark 2.3 and Lemma 2.4 we can easily find an N such that C m is invertible 
for all m > N , which implies that f m > 0. In this case by the properties of circulant 
matrices [0 the entries {/3k}™=- m+ i of C^ 1 can be computed as 



m— 1 



1 fm( 2m -l ) 



a = — E — 

A 2m- 1 f„+^ 

/— — m 

Now consider 

1/2 



= 27rifci/(2m-l) 



k = -m+l,...,m-l. (2.19) 



V 2 

' m-1 

/ _!_ e 27rife:E da; - ^ — e 27r 

/ /2 ^ l=-m+l /™(2^t) 



ii!/(2m-l) 



< 



1 



-1/2 



e 2 ™ fe;E dx 



1 ^ 1 

2m -1 ^ /(^-t) 6 

1 J V 2m — 1 / 



2-7riiW/(2ro-l) 



/ — — 771 + 



1 I 1 

— — y — - 

2m- 1 1 ^ ffs-i 



/ — — 771+ 



1 film-\) 



m— 1 1 

3 27T7fci/(2m-l) _ £_ 

f I 



3 2irifei/(2m-l) 



£=—771+ 



1 2 777-1 ) 



(2.20) 
(2.21) 



We estimate the expression above in two steps: 

1. We first consider ( |2~2(i| ). Note that = Efel-oo a k e 2 ™ kx and 



1 



E 



777— 1 / OO 



a fc e 



27r7fci/(2m-l) 



/( 2777-1 ) 

By setting S p = Y^L-oo a P +(2 m -i) q we get 



E ( E a P+(2m-i) 9 1 < 

p— — m+l \<J— — oo 



27T7ip/(27T7-l) 



^2m-l) p=-m+l 



8 p e 2 ^' {2m ^\ 1 = -m+l,..., m-1, 



and 



Hence 



1 m-1 1 
2m- 1 ^ . f(~£ 



5 -27rifep/(2m-l) 



+ 1 J V2m-P 



1/2 



e 2 ™ kx dx 



i m_1 1 



^lizikl I '(2m- 1) 



-1/2 



f(x) 2m -1 ^ f(Tr X -r) 

J v ' l=-m+l 1 v 2m- V 



| a* - < E l a fc+(2m-l)(/|- 
<t+0 



(2.22) 
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By Theorem |l.3|(a) A £ £ 1 implies A^ 1 S £ 71 . Hence we get for k = — m+ 1. . . . , m— 1 



l a fc+(2m-l) ? l 

9^0 



< 2ci e 

9=1 



- 7 i(fc+(2m-l)g) 



< 2ci e 
q=l 



7i((2m-l)g-m+l) 



2cie" 



7im 



1 _ e -7i(2m-l) ' 



(2.23) 



2. Now we estimate (2.21): 



m— 1 



1 I 1 

1 v - 

2m- 1 1 ^ ff^- 



= 27rifcZ/(2m-l) _ 



m— 1 

E 



27rifcZ/(2m-l) 



Z=-m+l 



i=_ m +l /(2m-l) 

1 '^C , '[ 



/m( 2m _i ) 



m— 1 



1 



2m 



rn — 1 



< 



2m 



1 , 1 

< max 



2m - 1 jir<m' 1 1 1 i= 5 +1 , £ 

< H^llllOl £ H- 

|fc|>m 



Z=-m+l /(2ra-l) /™(2m-l) 



(2.24) 



By Lemma 2.4 we can easily find for any e > an AT such that for all m > N there 
holds HC" 1 !! < (l+e)!^" 1 !!. Thus 



H^llllOl £ \a k \<(l + 8)1^-^2^1^ 

|fc|>m k=m 
p -im 

< ( l +£ )||A- 1 || 2 2 T -— - < {l + e)\\A-^fc 2 e-^ (2.25) 



for some 71 < 7 and some constant ci- By combining (2.22), (2.23), and (2.25) and 



hiding expressions as (1 + e)\A 1 1| in the constant c, we obtain estimate (2.1 



(ii): The proof of ( 2.1 7| ) is similar to the proof of (2.16). The only steps that 
require a modification are (2.25) and (2.25). By Theorem 1.3(a) A e Q s implies 
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A 1 e Q s . Hence we can estimate \ a k+(2m-i)q\ as follows. 

K+(2™-i)gl < c' 5^(1 + |(2m - 1) 9 + 

g#0 q^O 

OO 

= 2c'(l + 1 2m - 1 + fc|r s + 2c'^](l + | (2m - l)g + fc|)~ s 

g=2 

OO 

< 2c'(l + |2m- 1 -m + l|)~ s + 2c'^(l + |(2m- l)?-m+l| 

q=2 

OO 

< 2c'(l + m)~ s + 2c' / (1 + | (2m - l)x - m + l\)~ s dx 



s . (l + m)- s 2c(l + m)- s 

2c' l + m-s+20 ^ / < V ^ / . 2.26 
(2m — l)(s — 1) s — 1 



By adapting (2.25) to the case of polynomial decay we can estimate (2.21) by 



2 C ||A- 1 ||||0|^. (2.27) 



Combining (2.26) with (2.27) yields the desired result. □ 



The main portion of the proof of statement (ii) above is due to Gabriele Steidl [2J] . 

3. Error estimates for approximate solution of deconvolution problems. 

Consider the convolution of two sequences a = {afcj^L.oojC = {cfc}^L_ oc € ^ 2 (Z), 
given by b = a * c. Here a may represent an impulse response or a blurring function. 
Given a and b our goal is to compute c. This is known as deconvolution. In matrix 
notation the problem can be expressed as Lc = b where L is a biinfinite Toeplitz 
matrix with entries Lki — ctk-l- 

Sometimes a and c have compact support, in which case b also has compact (al- 
though larger) support and c can be computed by solving a finite banded Toeplitz 
system, see |M. It is well-known that this can be done efficiently via FFT by em- 
bedding the Toeplitz matrix into a circulant matrix. Of course, this approach is very 
attractive from a numerical viewpoint, at least if the system is well-conditioned (we 
discuss the ill-conditioned case in section ||) . 

However, if either a or c does not have compact support, the reduction of Lc = b 
to a finite Toeplitz system obviously will introduce a truncation error and embedding 
the Toeplitz matrix into a circulant yields an additional (perturbation) error. We 
nevertheless can try to make use of the FFT-based approach with the hope to get a 
good approximation to the solution. According to Q it has been shown in Q that 
the solution of doubly infinite convolution systems can be approximated by solutions 
of finite circulant systems. (Note that the approximation by finite circulant systems 
does not apply to one-sided infinite convolution equations.) 

Hence we are concerned with the problem of how good the approximation is 
obtained in that way and how fast the approximation converges to the true solution. 
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To answer these questions we proceed as follows. In the first step we study 
the approximate solution of a finite Toeplitz system by circulant embedding. In the 



second step we combine the obtained results with Theorem 2.1 to derive estimates for 
the approximate solution of doubly infinite convolution equations by finite circulant 
systems. The results that we will collect in the first step will also be very useful in 
section ^ for the analysis of circulant preconditioners. 

Let A n x^ = be given, where A n is an hpd nx n Toeplitz matrix. Of course, 
we have in mind that A n and y^ n > are finite sections of the biinfinite Toeplitz matrix 
L and the right-hand side y, respectively. 

As usual, we embed A n into a circulant matrix C^ n of size 2n x 2n as follows: 



C 2 



(3.1) 



where B n is the n x n Toeplitz matrix with first row given by bj- — a n -k for k = 
1, . . . ,n — 1. If a n is known we set bo — a n , otherwise we define bo = 0. For the 
following considerations it does not matter if we embed A n into a circulant matrix of 
size (2n — 1) x (2n — 1) or of size 2n x 2n (or larger). Choosing B n to be of the same 
size as A n is just more convenient for the proofs below. Since Ci n is circulant we can 
find A n again in the center of Ci n (and at any other position along the main diagonal 
of Cm), and express the embedding of A n into G% n as follows 





X 


X 


X 


Cm — 


X 


A n 


X 




X 


X 


X 



,(«) 



In spite of the biinifite system Lx — y in the background it is useful to embed 
symmetrically into a vector y^ of length 2n as follows 



y {n) = [0, 



0,2/"\0,...,0]. 



We assume for the moment that Ci n is invertible (and will justify this assumption 
later). An approximate solution z^ n ' to A n x' n > — j/W is now obtained by solving 
Cmz^ — y^ and setting — {i^™* 1 }^™ n %+i (^ e -' we ^ a ^ e as approximation the 
central part of z^ n > corresponding to the embedding of y^ into y^). 
We partition C^ 1 in the same way as Cin as follows 



°2n 



1 71 



where M n is a Toeplitz matrix and by Cauchy's interlace theorem (cf. Jl5[) invertible. 
Define S n := M^ 1 , then it is easy to see that z^ can be obtained as the solution of 

S n z {n) = y {n) . 

The question is now how well does z^ approximate iW. 
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Theorem 3.1. Let A n x^ — g/W be given where A n is annxn hermitian positive 

i 

2 

definite Toeplitz matrix with (A n )ki — a>k-i and au — J f(cu)e 27T ' lujk dco. Assume that 

_ i 

2 

\a k \ < ce-~*W and \(y (n) ) k \ < c'e-^ n/2 ~ kl . Suppose that C 2n as defined in (O) is 



invertible and let z^ n ' be the solution of S n z^ — y^ n \ where S n is the n x n leading 
principal submatrix of C^ 1 . Then there exists a 71 with < 71 < 7 and a constant 
C\ depending only on f mln and / max and on 71 such that 

\\ x (n) _^(«)|| < Cie -7in (3 3) 



Proof. Similar to Theorem 2.1 



we write 

\\ x W- z (n)\\<\\ A -i\\\\(A n -S n )S- l yM\\. (3.3) 

Using the Schur complement we can write S n as 

S n = A n -B n A- 1 B* n . (3.4) 

We will first show that the entries of the matrix A n — S n are exponentially decaying 
off the corners of the matrix. Note that equation ([TJ) implies 

A n — S n = E n 

where E n := BnA^ 1 B* n . We analyze the decay behavior of E n in two steps by 
considering first B n A~ x and then (B n A~ 1 )B^ l . 

Recall that \(A n ) k i\ < ce"^*-*! and note that \(B*)m\ < ce-^-^-^. By 



Theorem 1.3 we know that KA" 1 )^! < C2e _72 ' fc ~ i ' for some < 72 < 7, where C2 
depends on 72 and on / m ; n and / max , but is independent of n. We set S = 7 — 72. 
There holds 

n—l n—1 

|(i?^ 1 ) fe d = |^(B„) fc ,(4; 1 ) j; |<c C2 ^e"^"-l fc -^e-^^- i l. (3.5) 
3=0 j=0 

For simplicity we will absorb any constants arising throughout this proof that depend 



solely on 7 (or 72) in the constant c. We analyze the sum (3J5) further by splitting it 
up into three parts and in addition consider first the entries (BnA^ 1 )^ with k > I. 

(i) < j < I: 

l-l l-l 
c e -7(«-|fe-j|)g-72|j-(| _ ce ~f («-fc)g-72( g-(7-72)j < cg -72(n-fc+0 

3=0 3=0 

(ii) l<j<k: 

-7(n— |fc— 3\) <r ^ c -l( n ~ k ) 

1 — e -(7+72) 



E e -7(«-|fc-il)g-72b'-'| < C g-7(™-fc) e 72i_r < ce -72(n-fc+0 

— 1 — p-(7+72) — 



3=1 
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(iii) k < j < n: 

11-1 P (7-72)(n-l) 
c ^ e -7(n-|fc-j|) e -Ti|j-'l < ce -T(™+fe) e 72i _ (7 _ 72) (3.6) 

„5(n-l) 

= ce -72(™+fc-0 e -<5(n+fc)_r < ce -72(«-i+fc) < ce -72(n-fe+0 

1 — e -(7~72) — — 

Similar expressions can be obtained for the case I > k by interchanging the roles of k 
and I in the derivations above. Thus 

We now estimate the decay of the entries of B n A~ x B n . Since B n A^ x B^ is her- 
mitian, it is sufficient to consider only the entries (B n A~ x B^)ki with k > I. There 
holds 

n-l 

\{B n A- x B* n ) kl \ < C ]T e -72(«-l^l) e -7(«-b-*|). (3 . 7) 

3=0 

As before we proceed by splitting up this sum into three parts, 
(i) < j < I: 



i-i 

c \~^ e -J2(n-\k-j\) e -"/(n-\j-l\) _ ce -(7+72)" e 72fc e 7' < C g-72(2n-fe-() 

^ 1 — e -(7+72) — 

3=0 



(ii) l<j<k: 

e -l2{n-\k-j\) e -y{n-\j-l\) < ce -( l2+ 5)(n+l) g- 7 2(n-fc) £ < ce - 72 (2n-fc+;) 

- 1 _ e -5 - 

(iii) k < j < n: 

n ~ X =(7+72)(n-l) _ ( 7+72 )(fc_l) 



C 



e -72(n-|fc-j|) e - 7 (n-|j-t|) < g-(7+72)« e -72(fc+0 f ! 

j=k 

< ce -72(fe+0. 

Hence, by combining (i), (ii), and (iii) we get 

\(B n A~ x B*) kl \ < c (e- 2 ^ n e^ k+l) + e - 2 ^ n e^ k -^ + e~^ k+l ^). (3.8) 
The entries of Ci n satisfy 

,, x , \ce-^ k - l \ for =0,...,n- 1, , 

GhOifc < < , „ m (3-9) 
ce - 7 (2n-|fe-i|) for |fc - Z| = n, . . . , 2n - 1. 
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By Theorem 2J3 there exists a 73 < 7 and a constant C3 depending on 73 and on f n 
and / max such that 



2n Ilk 



< 



c 3 e 
c 3 e 



-1s\k-l\ 



for \k - l\ 



- j3 (2n-\k-l\) for | fe _ 2 | 



0. 



n, . 



,n - 1, 
,2n-l. 



(3.10) 



Hence z^ := S~ l y^ satisfies < c 3 e~^ n / 2 ~ k ^. Set u< n ) = (A n -S n )z^ n l After 
some lengthy but straightforward computations we get 



\(u (n) ) k \ < 



Cie - 7l (n/2+fe) for fe = 0,...,n/2, 
Cie -7i(3n/2-fe) for k = n/2,...,n 



(3-11) 



for some 71 < 73. Hence ||w^ n - ) |j < ce -73 ™/ 2 , which together with (3^3) completes the 
proof. □ 

Corollary 3.2. Let Lx — y be given where L is a biinfinite hermitian positive 



definite Toeplitz matrix with entries Lki = ak-i and let S n be as defined in (3.4). 

Assume Le£ 7 and \yk\ < e _7 ' fc ' and let z^ be the solution of S n z^ = . Then 
there exists an N such that for all n > N 

Wx-z^W < ce-^ in . 



for some < 71 < 7 and a constant c independent of n. 

Proof. First note that by Remark [2.3] we can always find an TV such that S n exists 
for all n > N. There holds 

\\x - = WL-'y S-y«)|| < WL-'y - L?y™\\ + \\L?yW - S^y^W (3.12) 

where A n is an n x n finite section of L. The result follows now by applying Theo- 



rem 2.1 and Theorem 3.1. □ 



Theorem 2.1 



and Corollary 3.2 provide two different ways to approximate the 



solution of biinfinite Toeplitz systems. Which of the two is preferable? This depends 
on two criteria: (i) The accuracy of the approximation for given dimension n; (ii) the 
computational costs for solving each of the finite-dimensional systems. 



The solution of the circulant system in Corollary 3.2 can be computed via 3 



FFTs of size 2n. The Toeplitz system in Theorem |2.l| can be solved by the conjugate 
gradient method in approximately 3/c FFTs of size 2n, where fc is the number of 
iterations. Of course, additional preconditioning can significantly reduce this number 
at the cost of two additional FFTs per iteration (see also Section ||). For both, the 
circulant and the Toeplitz system, zeropadding can be used to extend the vectors to 
"power-of-two" -length. 

Example 2: We consider the same biinfinite Toeplitz system as in Example 1. 
We compare the error when approximating the solution by using the circulant system 
of Corollary [O and by the Toeplitz system of Theorem 2.1. We compute for each 
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10' 



error via Toeplitz system 

- - error via circulant system 




150 200 
dimension n 



350 



Fig. 3.1. Comparison of error for the solution of a biinfinite Toeplitz system with poly nomial 
decay. We compare the approximati on e rror of the Toeplitz system described in Theorem 
of the circulant system of Corollary 



2.1 



to that 



for increasing matrix dimension. The approximation error 
of both methods is almost identical, so that the difference between the two graphs is hardly visible. 



1, 350 the approximation error || x — x^ || and ||x — z^ n '\\ respectively. As can 



be seen from Figure 3.1 both methods give almost the same error, in fact the two lines 



showing the error are hardly distinguishable. A similar behavior can be observed for 
other examples involving biinfinite Toeplitz matrices with fast decay. Since solving a 
circulant system is cheaper than solving a Toeplitz system, the approximation scheme 



of Corollary 3.2 seems to be preferable in such situations. 

Many variations of the theme are possible. For instance if A n is an s-banded 
(biinfinite) Toeplitz matrix with s < n/2, we could use Strang's preconditioner as 



approximate inverse. Due to the explicit constants in Theorem 1.1 this approach 
allows us to give an error estimate with explicit constants (cf. also Theorem 5 in p5i ) . 



Theorem 3.3. Let A n x^ — y( n > be given where A n is an n x n hermitian s- 
banded Toeplitz matrix with s < n/3 and positive generating function and let = 
for | n/2 — k\ > s. Let C n be the n x n circulant matrix with first row given by 
(an, ffli, . . . , a s , 0, . . . , 0, a s , . . . , ai) and let z^ n ' be the solution of C n z^ = y( n > . Then 

\\ X W - < 3V2c\" in (\-' 1s - A" 7 ( s+1 ))- 3 



where c and A are as in Theorem 1.1 



Proof. The proof is similar to that of Theorem [O . To avoid unnecessary repeti- 
tions we only indicate the modifications, that are required. 
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By Remark 2.3 C n is invertible. Note that C n is a matrix with three bands, one 
band is centered at the main diagonal, and the two other bands of width 2s are located 
at the lower left and upper right corner of the matrix. It follows from Proposition 5.1 
in H that the entries of C~ x decay exponentially off the diagonal and off the lower 
right and upper left corner. More precisely, 



I(OmI< 



c\\ k - l \ ifO<|A;-Z|<n 
cA 2n+i-|*-J| if n+ i< \k-i\ < 2n. 



where c is as in Theorem 1.1 with A = . With this result at hand it is easy to show 
that the entries of C~ 1 y^ decay exponentially. 

A n — C n has a simple form, it is a Toeplitz matrix with first row 

(0, . . .,0,a g o a _i, . . . ,ai). 

When we compute vS n > := (A n — C n )(C~ 1 y( n * > ) the non-zero entries of (A n — 
C n ) are multiplied by exponentially decaying entries due to the exponential decay of 
((C-y n) ). This leads to the estimate ||u< n )|| < 2^/2c{\ s - A S+1 )- 3 A". □ 

An interesting alternative to periodic boundary conditions is the use of Neumann 
boundary conditions considered in ppf . This modification will be discussed elsewhere. 

4. Preconditioning by embedding and exponentially decaying Toeplitz 
matrices. The accuracy of the solution of the Toeplitz system A n x^ — y( n ' by 
using S~ 1 as approximate inverse depends crucially on the decay properties of the 
right hand side yW, If y( n > does not have appropriate decay conditions the approach 
in section ^ may not yield an approximation with sufficient accuracy. But we can still 
use S 1 " 1 as preconditioner and solve by the preconditioned conjugate 

gradient method 0] . The construction of preconditioners via circulant embedding is 
well known, it has been thoroughly investigated in |6j and for the special case of band 
Toeplitz matrices in Q . 

Compared to the certainly more elegant and more general approach in ||, the 
approach undertaken in this section has the advantage that it yields some quantitative 
results. It shows that the clustering behavior of the preconditioned matrix S~ 1 A n is 
the stronger the faster the decay of A n is. Moreover, our approach allows us to prove 
a conjecture by Nagy et al., cf. [fill], and will provide a theoretical explanation for 



some numerical results presented in [19[ ] and 14 . 

The theoretical analysis of the clustering behavior of the eigenvalues of the pre- 
conditioned matrix S~ 1 A n is inspired by the work of Raymond Chan |S|. We will 
show that S~ 1 A n can be written as S^An = I n + R n + K n , where /„ is the identity 
matrix, R n is a matrix of small rank, and K n is a matrix of small 2-norm. 

Theorem 4.1. Let A n be a hermitian Toeplitz matrix whose entries decay 
exponentially, i.e., |dfc| < Ce 7 ' fc ', k = 0, . . . ,n — 1. Set S n = A n — BnA^ 1 B* t , where 



B n is as defined in (3.1). Then for all e > 0, there exist N and M such that for all 



n > N at most M eigenvalues of A n — S n have absolute value exceeding e. 
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Proof. By definition of S n we have 

A n S n — E n 



where E n = BnA^ 1 B* n . We know from equation ( |3.8| ) in the proof of Theorem 3.1 
that the entries of E n can be bounded by 

|(JS^)kj| < c ( e - 2 Ti"e 7l(fc+i) +e- 27l "e 7l|fc ~' 1 +e-' 1l(k+l) ). (4.1) 

For x = [xo, Xi, . . . , x n -2, x n -i] G C" and N < n/2 we define the orthogonal 
projection P^ by 

Pnx = [0, ... ,0,x N ,x N+ i, . . . ,x n _jv-i,0, ... ,0], 

and identify the image of P N with C n - 2N . We set E n N) = P N E n P N . In words, E { n N) 
is obtained from E n by taking only the central (n — 27V) x (n — 2N) submatrix of E n 
and setting the other entries surrounding this block equal to zero. Then E n - E n N) 
has 2N "full" rows and n — 2N "sparse" rows, where each of the latter rows has 
non-zero entries only at the first N and last N coordinates. Thus the dimension of 
the space spanned by the sparse rows is at most 2N. Hence rank(£' n — En" 1 ) < AN. 
Due to the decay properties of E n it is easy to see that 

n-N-l 



Nl I 
l=N 

Using ( |4.lD we get after some straightforward calculations 

\\E n N) ||i < c(e~ 71 " + e -^ n - 2N) + e- 2 ^ N ). (4.2) 

It is obvious that for each given e > we can find an N such that for all n > 2N 
H-En^lli < £■ Since E n N ^ is hermitian, we have \\E^ N '\\i = WE^ W^. Thus 

ll^lb < (H^Wlllll^lloo) 1 / 2 < £ - 

Hence for large n the spectrum of En lies in (— e,e). By the Cauchy interlace 
theorem we conclude that at most AN eigenvalues of A n — S n have absolute value 
exceeding e. □ 

Lemma 2.4 implies that for any s > we can find an M such that for all n > M 
S n and S" 1 exist and US'" 1 !! is bounded by / m i n — £ > 0. Proceeding as in ||, Chapter 
2, we express 5~ 1 A n as 

1 A n = I n + S n 1 (A n — S n ) 

and arrive at the following 

Corollary 4.2. Let A = {a^} be a hermitian positive definite Toeplitz matrix 
with |dfe| < ce" 7 '' ' for 7 > 0. Then for all e > there exist N and M > suc/i i/iai 
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for all n > M at most N eigenvalues of S n A n — I n have absolute values larger than 



With the results presented in this paper it should not be difficult for the reader 



to derive Theorem 4.1 and Corollary 4.2 for Toeplitz matrices with polynomial decay. 



Theorem 4.1 (in particular fl4.2} )) and Corollary 4.2 show that the clustering behavior 
of S~ 1 A n is the stronger the faster the decay of the entries of the Toeplitz matrix is. 

In Jll| and in jTJj Nagy et al. consider the solution of convolution equations and 
banded Toeplitz systems by the preconditioned conjugate gradient method using a 
preconditioner similar to the one in this section. In ]l9fi the Toeplitz matrix is rectan- 
gular, but the embedding is in principle the same. The theoretical results presented 
in []l9[ [l4| only hold for banded Toeplitz systems, but in the numerical experiments 
the authors consider also non-banded Toeplitz systems, where the Toeplitz matrix has 
exponential decay (see jl9|]) or polynomial decay (see Example 4 in juj). It is noted 
in jl9| that "it is surprising that the number of iterations is still quite small". The 
authors also point out that the numerical experiments indicate additional clustering 
around one of the spectrum of the preconditioner matrix, which is not covered by 



their theoretical results. With Theorem |hlj and Corollary [4.2| at hand we can provide 
a theoretical explanation for the numerical observations in Jl!| and |Q, at least for 
the 1-D case. The key lies in the fast decay of the inverse of the Toeplitz matrix. 
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